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ABSTRACT 



Supersonic flow past harmonically vibrating two- 
dimensional panels and cylindrical shells is analyzed 
using a linearized method-of-characteristics procedure. 
A detailed description of this method to solve the 
linearized unsteady potential equation is given and 
numerical results are presented to indicate the nature 
of the aerodynamic pressure distributions. Also, com- 
parisons of the present results are made with earlier 
work by Nelson and Cunningham and by Anderson which is 
based upon quite different approaches. 
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I. 



INTRODUCTION 



The determination of the pressure distribution over 
vibrating panels and shells is a prerequisite for the study 
of aeroelastic stability. This problem has been previously 
considered by a number of investigators, who applied either 
various asymptotic approximation techniques (e.g., piston 
theory) or operational methods to the solution of the lin- 
earized unsteady potential equation. For a recent compre- 
hensive review of this subject refer to Reference 7. In 
the present study a quite different approach is developed; 
i.e., the method of characteristics is used to obtain a 
solution to the problem of supersonic flow past vibrating 
panels and shells. 

The major assumptions in the mathematical model are 
that we are dealing with a perfect gas, that the flow is 
irrotational and supersonic, and that linearization is 
permissible. Reference 1 was used to help develop the 
linearized unsteady potential equation and the boundary 
conditions for panels and shells. These equations were 
then rewritten in cylindrical coordinates . 

Steady supersonic flow past bodies of revolution was 
studied by several authors using linearized characteristics 
methods; i.e., Haack [Ref* 9], Sauer [Ref. 1 6 ] , Oswatitsch 
and Erdmann [Ref. l£J . Teipel [Ref. 12} developed a char- 
acteristics procedure for two-dimensional supersonic flow 
past airfoils oscillating at arbitrary frequency using the 
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velocity perturbations and the velocity of sound perturbation 



extension of the Oswatitsch-Erdmann method to slowly oscil- 
lating bodies of revolution. 

In the present report a linearized characteristics method 
is presented to study supersonic flow past cylindrical shells 
vibrating at arbitrary frequency. The velocity potential 
function is introduced as the dependent variable and the 
linearized unsteady potential equation is solved subject to 
the proper boundary conditions. Although two-dimensional 
flow past vibrating panels is contained in this solution as 
a special case for infinite shell radius, two separate 
programs were written in Fortran IV code for the IBM 360 
system, one for panels and one for shells. Thus, the panel 
program is complementary to Teipel ' s work in that it uses 
the velocity potential rather than the velocity perturbations 



served as an independent method to verify the results. 
Anderson's work, on the other hand, was used to compare the 
method-of-characteristics predictions for cylindrical shells. 



as dependent variables. Platzer and Sherer 




gave an 



as dependent variables. The method of singularities 
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II. PROBLEM FORMULATION 



In nearly all engineering studies some basic assumptions 
have to be made to formulate the problem. In this study 
this assumption is that the compressible medium with which 
we are dealing is a perfect gas. We therefore have the 
equation of state for a gas 

P = pRT (2.1) 

Flow of a perfect gas is then completely described by 

% 

expressing the ambient quantities P,p and T as well as the 
velocity components u, v, w as a function of position x, y, 
z and time t. These six dependent variables are always 
related by the conditions of conservation of mass, momentum 
and energy as developed in Johns [Ref. l] . 

Conservation of Mass 

9 (pu) 3 (pv) 3 (pw) 3p _ 

3x + 3y 3z at 

Conservation of Momentum 

1 3P 3u , — 3u , ,, 3u , TI 3u 
p 3x d t 9x 8y 9z 

1 3P 3v _ 3v 3v 3v 

" p 37 - 3t + U 3^ +v 3y +w 3^ 

1 3P 3w — 3w 3w 3w 

- — T— = -r-r + U -T— + V -T— +W TT 

P dZ d t dX dy d Z 



( 2 . 2 ) 



(2.3) 



/ 



10 



Conservation of Energy 



§t f 5 



W 



1 

P 



g g 

3x (^ p ) + ay” ( vE> ) + 



7£( wp )] 



where 



W* = + v 2 + w 2 . 



1 2 = u 2 

Equation 2.4 can also be stated as (See Ref. 1) 



(2.4) 



(2.5) 



PS 

Dt 



0 



( 2 . 6 ) 



The condition of conservation of energy therefore reduces 
to the simple statement that the entropy of a fluid particle 
remains constant. For later convenience we now rewrite the 
continuity equation, Eq. 2.2 



Dp 

Dt 



+ 





0 



Introducing the speed of sound 



(2.7) 



a 



2 




( 2 . 8 ) 



we have 

Dp _ 9 p DP _ 1 DP 
Dt 3P Dt 2 Dt 

Cl 

Therefore the continuity equation can also be written 



DP _ _ D a 2 ( 9u 8v 3w\ 
Dt [JZ + 97 + TiJ 



(2.9) 



The velocity potential is introduced to reduce the 
number of dependent variables. Its existence depends upon 
the condition of irrotationality of the flow. This condition 
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is expressed mathematically by the vanishing of the curl 
of the velocity vector and can be written in component form 



civ jha 9w 9>v _ 9_u 

3x “ 9y = 0' 9y ” 3z = ^ • 9z 



3w 

3x = °* 



( 2 . 10 ) 



Ref 



•’] 



that flows 



It can be shown by Kelvin's theorem 
starting either from reservoirs or from parallel stream- 
lines and being irrotational initially will remain so. Only 
strong shock or intense heat will require a reformulation 
of this law. However, such effects are usually quite 
unimportant in aeroelastic considerations. Therefore, Eq . 
2.10 is necessary and sufficient to derive the velocity 
components u, v, w from the velocity potential cf> such that 



9<J> 

u = 9x ' 



3 (J) 

v " 3y 



9(f) 

W = -sr-*- 

9z 



( 2 . 11 ) 



Therefore the continuity equation can be written 



1_ Dp _ 9u 
p Dt ~ 3x 



+ _ + ^ - lit + i!± + 



9v 

3y 



9w 
3 z 



2 

3z 



2 



3x 3y 

In vector form the Eulerian equations of motion can also 
be written as 



( 2 . 12 ) 



3v 1 

Yt + (v * grad)v = - ± grad P. 

This is equivalent in irrotational barotropic flow to 



(2.13) 



at 9 rad $ + 9 rad \ y ~ = 



w 



- — grad P 
P 



(2.14) 
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or 



grad 



3 <J) W 2 [dp] 

[dt + — + Jp“J= 0 



Integration therefore gives 

r 



3<f> W" 

te + r 



dP 

7T- = C(t) 



(2.15) 



J 



Moreover, C(t) = i U 2 whenever the remote fluid motion 

2 

consists of parallel streamlines with velocity U. The 
partial derivative of Eq. 2.15 with respect to time is 



3 r 3 4> W 2 ] 3 

3 1 (_3t + Tj ~ 3t 






dP 

P 



1 3P 
p 3 1 



(2.16) 



In rewriting the continuity equation 



d> + <4,7,7 "f 4* 
Y xx yy zz 



(2.17) 

a z p 



one obtains therefore as a result of Eqs. 2.3 and 2.16 



1 - 



x 



, 4 + 

,2 l V XX 



1 -- 



*TJ ^-TT 



6 + 
r zz 



2( My 



26 6 
y x Y y 



2 4> y Z 



r xy 



'xz 



*yz + 



(2.18) 



24 



x 



,2 Y X t 



24, 



24. 



. 2 Y yt 



J Vzt " ^tt 



= o 



Let a denote the speed of sound. For infinitely large 
speeds of sound (incompressible flow) the above equation 
reduces to Laplace's equation. Equation 2.18 is the 
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governing differential equation for general nonstationary, 
nonviscous, potential flow and is not limited to small 
disturbances. However, because of its nonlinear nature, 
solutions were found only in very few special cases. There- 
fore, it is necessary to introduce the small disturbance 
concept. Regarding all velocity disturbances small in 
comparison with U, a, and (U - a) and pressure differences 
and density changes small in comparison with main stream 
pressure and density, a perturbation potential is introduced 



If a is expanded around its value for the undisturbed stream 
the procedure of retaining only first order terms leads to 
the following equation for the disturbance velocity potential 



where c is the constant value of the speed of sound in the 
uniform stream and M the constant free-stream Mach number. 

It is this linearized unsteady potential Eq. 2.20 which is 
being used in most aeroelastic analyses as the basic equation 
to predict the oscillatory pressures and, generalized aero- 
dynamic forces. However, it must be emphasized that its 
validity is restricted to purely subsonic or supersonic 
flows, as well as to sufficiently unsteady transonic flows. 



so that 



<(> = Ux + 



(2.19) 




( 2 . 20 ) 
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Expressing Eq. 2.20 in cylindrical coordinates one has 






XX 



c? 



rr 



+ r ^ + ^7 ?° e ' jh [^tt + 2u fet + 



= 0 



( 2 . 21 ) 



There are two main boundary conditions applicable to a 
vibrating cylinder in supersonic flow. The first is that 
the panel surface is impenetrable to the gas, i.e., the flow 
is tangential to the surface at each instant of time. The 
second condition is prescribed by the Sommerfeld radiation 
condition which requires waves to propagate away from sources 
of disturbance toward infinity. Furthermore, uniform super- 
sonic flow must exist upstream of the panel leading edge, 
since no disturbances can propagate upstream. 

If the equation of the surface of a body which is moving 
in a time-dependent manner is given by 



F(x,y,z,t) = 0 (2.22) 

then it can be shown that the flow tangency condition 
requires [Ref. 8) 



DF 9F 
Dt “ 8t 



- 

9x 



9F 

9y 



+ xi ~ + v — + w — = 0 



9F 
9 z 



(2.23) 



This condition states that the rate of change of the numerical 
value of the function F is zero, when we follow the motion 
of a particular fluid element (substantial derivative) , so 
that the element continually touches the surface F = 0. 
Consider a two-dimensional panel (Fig. la) whose chord is 
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aligned with the x-axis and whose upper surface is exposed 
to a uniform flow of velocity U in the positive x-direction. 
We then find that we can describe the body by 



F = z-h(x,t) = 0 



(2.24) 



. 3 F 

Since — = 1 we have from the boundary condition 
3 z 



w _ 3h — 3h 

w ' at + u 35? 



For z = h(x,t) 



(2.25) 



Imposing the linear perturbation assumption Eq. 2.19 we have 

<t> = Ux + <p (2.19) 

where the disturbance velocity components are obtained from 



u = 91 



w 



_ 3<P 



3x 3 z 

and satisfy the order-of-magnitude requirement 



(2.26) 



u, w << U 



(2.27) 



dh 



The terms u — can then be neglected compared to the much 



r\ 

larger term U — . This leads to 



3h 



3h 



W ' 3 1 + U 3x 

It can be shown that it is possible to express the flow 
tangency condition in the plane z = 0. Taylor series 
expansion gives 

w(x,z,t) = w ( z , 0 + , t) + h 8W(X,0+,t) 



(2.28) 



+ , 



h 2 3?w(x,0 + ,t) 



2« 



3z 



+ ---- + 



(2.29) 



3 Z‘ 
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X 



FIGURE la TWO-DIMENSIONAL PANEL 




FIGURE lb PROFILE AND CROSS SECTION 
OF CYLINDRICAL SHELL 
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Using the same argument as above the product terms are 
neglected so that the flow tangency condition equals 



w 



8h 

3t 



+ U 111 
3x 



for z = 0 + 



A similar consideration for the cylindrical shell 
whose vibration mode is given by 

iuit 



h(x,e,t> = z < x > Cos ne e 

results in the following boundary conditions 



(Fig. 



3 r 



r=R 



U 





Cos n© 



iwt 

e 



= 0 



for 0 < x < L 
for x < 0 



The calculation of perturbation pressures and pressure 
coefficients begins with Eq. 2.15 




For uniform parallel flow far upstream we have 
C(t) = 2 . 

Using again the perturbation assumption, Eq. 2.19 

<J> = Ux +<p 



and noting that 

W 2 = |u + 
we can approximate 



dj£ 

ax 




H + 

at 



w 2 

2 



d£ 

at 



+ 



Uj 

2 



+ U 



a_^> 

ax 



(2.30) 
lb) 

(2.31) 



(2.32) 



(2.5) 

, ' 

(2.33) 
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by dropping second order' terms. Also the density-pressure 
relation is expanded as 




(2.35) 



Hence, the perturbation pressure is related to the perturbation 
potential by 



P-Pco = " Pco- 



3 3 

it + 0 H 



(2.36) 



Writing this equation in pressure-coefficient form one has 



Cp 




2_ 3<P _ 2 

U 2 at U 9X 



(2.37) 



Consistent with the assumed cylinder deflection, Eq. 2.31 a 
velocity potential is now assumed in the form 

= <3>(x,r) Cos n0 e^ ut (2.38) 

Substituting this into the Eq. 2.21 results in 

(1-M 2 ) $ xx + $ rr + ~ - pz a z $x + = ® (2.39) 

This is the basic equation for <Mx,r), which needs to be 
solved, subject to the boundary conditions 
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r=R 



For 0 < x < L 



ii 

9r 



= u 



az 

9x 



+ lccZ 



= 0 



For x < 0 



(2.40) 



yr 27 

When non-dimensional coordinates x = jr, r = £ and the 

(oL 

reduced frequency K = jj— are introduced then it is found 

(1-M 2 )$ +4 + 7 ® $- 2iKM2 * + K2M2 $ = 0 (2.41) 

xx rr £ 2 x 



The bar has been omitted for simplicity. The boundary 
condition is 



$ 

*r 



r=R 



tf + iKZ ] 



= 0 



0 < x < L 

x < 0 (2.42) 



Converting Cp to cylindrical coordinates and using the above 
non-dimensional coordinates one has 



Cp 



-2 iK$ 







Cos nOe 



(2.43) 
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III. NUMERICAL TECHNIQUE 



A. METHOD OF CHARACTERISTICS 

Consider the following pair of simultaneous first order 
partial differential equations for the dependent variables 
u and v. 



1 ax 



+ B. 



+ 

ay 



+ Ui 



ay 



(3.1) 



A ? + 

2 ax 



B ( 



3u 



9Z = 0 

ay 



where the coefficients A^ , . . . , are functions of u 

and v but not of x and y. 




FIGURE 2 
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Assume that in Figure 1 the solution is known up to the 
curve CPC. At P one knows continuously differentiable 
values of u and v along the curve CPC as well as all their 
derivatives in the directions that lie below the curve. 
However, the question needs to be answered whether there is 
enough information to deduce the directional derivatives 
of u and v at P in directions that lie above the curve. 

The directional derivative of u in the direction s is: 

8u _ 8_u 8_x 8_u 8_y 

8s “ 8x 8s + 8y 8s (3.2) 

so that the directional derivative in any direction is known 
as soon as the derivatives 8u/8x and 8u/8y are known. 

To find the values of 8u/8x and 8u/8y, write out 
equations 3.1 at point P and the directional differentials' 1 ' 
of u and v taken along CPC at P. In matrix form the four 
equations are: 



A 1 


B 1 


c i 


d T 




Tu/8x 




D " 


A 2 


b 2 


C 2 


°2 




8u/8y 




0 


d x 


. dy 


0 


0 




8v/8x 


= 


du 


0 


0 


dx 


dy_ 


1 


8v/8y 
— » -- 


. 


dv 



With u and v known at P the coefficients A 1 , . . . • > D 2 

will also be known. When the' directions of CPC are known 
then dx and dy can be found. If u and v are then known along 



1 In Eq. 3.3 du is an abbreviation for 8u/8s ds, where 
3 u/3 s is given by Figure 3.2 and s now measures distance _ 
along the curve CPC. The symbols dv , dx , and dy stand for 
similar abbreviations. 
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CPC, du and dv will be known. Equations 3.3 therefore 
constitute a pair of four simultaneous linear algebraic 
equations for the four first derivatives. If the determinant 
of this matrix is not zero, there is a unique solution. In- 
this case satisfaction of these governing equations requires 
that the directional derivatives have the same value above 
and below the line CPC. If the matrix determinant had 
equaled zero, there would not be a unique solution and 
there could be discontinuities across the line CPC. 

Expanding the determinant of Eqs. 3.3 and setting it 
equal to zero one has 

(A C^A^^dy 2 - (A 1 D 2 -A 2 D 1 +B 1 C 2 -B 2 C 1 ) dxdy + 

• (B 1 D 2 -B 2 D i )dx 2 = 0 (3.4) 

This is a quadratic equation for the slope dy/dx. If the 
direction of the curve CPC at P is such that it has a slope 
satisfying Eq. 3.4, then the derivatives of u and v are not 
uniquely determined by the values of u and v along the curve. 
Such a direction is called a characteristic direction. The 
quadratic Eq. 3.4 gives two real slopes, one real slope or 
a pair of complex slopes, depending upon whether the 
discriminant 

( A^D 2 "”A 2 D^tB ^C 2 "”B 2 C^ ) 2 — 4 (A^C 2 “A 2 C^) (B^D 2 “B 2 D^) (3 • 5) 

is positive, zero or negative. This is also the criterion 
for cataloging the system Eq. 3.1 as hyperbolic, parabolic, . 
of elliptic. The system is hyperbolic at a point if there 
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are two real characteristic directions; it is parabolic if 
there is only a single characteristic direction; or it is 
elliptic if there are no real characteristic directions at 
the point. 

The foregoing analysis can be repeated for the single- 
order quasi-linear equation. 



a 



9 ^ 



+ b 



ifi . 

9x9y 



c 



9 2 ^ _ , 
ay z “ f 



(3.6) 



in which a, b, c, and f are functions of x, y, ip, , and 
dip/dy. The characteristic directions are determined from 
the quadratic 



ady 2 - bdxdy + cdx 2 =0 (3.7) 

and hence Eq. 3.7 is hyperbolic, parabolic or elliptic at 
a point according to whether the discriminant b 2 - 4ac is 
positive, zero, or negative. 

To determine the characteristic curves we return to the 
pair of first-order Eqs. 3.1 and suppose that the system is 
hyperbolic throughout the domain of interest. At every point 
there are two roots, (dy/dx)a and (dy/dx)B , to the quadratic 
Eqs. 3.4. A curve, which at each of its points has the slope, 
(dy/dx)a , is said to be a a characteristic; A curve, whose 
slope is everywhere, (dx/dy)g, is said to be a $ characteristic. 
Thus there are two families of characteristic curves filling 
the domain as shown in Figure 3. 
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FIGURE 3 



If we are considering a characteristic direction so that 
the determinant in Eq. 3.3 is zero, the right-hand column 
must be compatible with this if there are to be any solutions 
at all for the first derivatives; i.e., when the right-hand 
column is substituted for any of the columns on the left, 
the resulting determinant must also vanish. Replacing the 
fourth column on the left with the column on the right and 
setting the determinant equal to zero leads to 



Insert (dy/dx)a into Eq. 3.8; "it then becomes an ordinary 
differential equation of u and v along the characteristic. 

A similar equation can be obtained along a 8 characteristic. 




dv = 0 . (3.8) 
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A method of attack has been outlined for solving hyper- 
bolic systems: first, locate the characteristic curves; 

second, integrate the ordinary differential Eq. 3.8 along 
the characteristic. This procedure is called the method 
of characteristics. 

B. FORMULATION OF THE EQUATIONS PROGRAMMED 

In Section II, D, Eq. 2.41 was developed as the partial 
differential equation for $(x,r) that needed to be solved 
subject to the boundary condition Eq. 2.42. If M > 1, the 
preceding equation is of the hyperbolic type and possesses 
real characteristics, which satisfy the ordinary differential 
equation 



(M 2 - 1) dr 2 - dx 2 = 0 



(3.9) 



Introduce along the characteristics the arc-length ds 2 =M 2 dr 2 

, , , , dx /M 2 — 1 i _ dr _ 1 /o i a \ 

then one has X = = — rr — , r - -r— - _ VJ.iu; 

ds M ds m 

Note: The angle between the characteristics and the x-axis 

is Sin a = — (3.11) 

M 

An arbitrary function F(x,r) has the following derivative 
along a characteristic 



dF 

31 = F x X ’ + F r 



r i _ A p 4- p 
r M x " m r 



(3.12) 



Denoting the derivative along the left-running characteristic 
with F^ , the derivative along the right-running characteristic 



with F 2 one has 
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(3.13) 






F + — F F = 

v M r r* ' r 2 



M x M r 



y/^-l 

M 



F - — F 
x M r 



Solving for F x , results in 



F r “ f < P 1 - F 2>' F = 



M 



(F, + F 9 ) 



X 2/m 2 -1 



(3.14) 



Call s^ the arc-length of the left-running characteristic, 
then one has for the second derivatives 



12 



dF 1 
dS^ ' 
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dF; 

dST 



(3.15) 



Hence 



12 



= (F-i ) 



/m*-1 



l'x 



M 



- (F l>r 



1 

M 



and 



-F 



12 



= i- [d- 

M 2 t 



M2)F X x +F rr| = - F 21 



0 



Therefore the basic equation for $(x,r) was written 



■12 ■ 57 l'l 2 



-M 2 $, „ + ($ n _ ihM 3 ($ 1+ $ 2 ] + (k 2 M 2 - jC) $ = 0 



/M 2 -! 



(3.16) 



or 



* 12 -i^r ( $1 + * 2 ) + ( K * ' ^ * 



This equation was then written in finite difference form. 

Using the backward difference form of the finite dif- 
ference method as developed in Wang Ref. 3 one obtained 



3x 



) j " dr +c Pj) Applying 



this technique to 
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Eq. 3.16 where the point's A and B were on the right-running 
characteristic (See Fig. 3) then 



^ (B) - ^ (A) 




4> (3.17) 



AS 




FIGURE 4 



The values of and in Eq. 3.17 were not specified at 



a point. When this equation was first written these values 
were taken at A, the previous point. It was found in 
comparing the results in Chapter V that there was a consider- 
able amount of error present. In investigating the source 
of these errors it was found that if the values for <i>^ and 
$2 were averaged between points A and B that these errors 
disappeared. The value for $ at point B was not yet known 
but it was found by integrating along the characteristic. 
Using trapezoidal rule 



Substituting this into Eq. 3.17 where $ (B) arises, and using 

, ^ 

the above averaging technique for $2 and $ results in 




(3.18) 
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A T, AS iKM ~1 T AS iKM 1 

^1 (B) L 1 " 4r (B) M + 2/M 2 -1 ASj =$ 1 (A) |l+ 4r (A)M " I / M 7 -! AS J " 



*2 (B) 


r as 


IKM 


AS 


- AS2 (k 2 


n 2 \ 






|_4r (B) M 


+ 2/M 2 -l 


4 [ K 


r 2 (B) M 2 J 






$ 2 (A) 


r as 


iKM 


AS 


AS 2 / 


h2 \~| 






[_ 4r (A) M 


+ 2 /M 2 - 1 


- — ^K 2 - 


r 2 (A) M 2 jj 


+ 


(3.19) 


$ (A) 


1*. - , 


n 2 1 AS 


t 


r- 

$ (A) K 2 - 


n 2 "j AS 




L r 7 


(A) M 2 J 2“ 


T 


r 2 (B)M 2 J 2 




Similarly to 


solve for 


$ 2 


at point B, 


the above 


procedure 



was used except now move along the left running character- 
istic and integrate along the left running characteristic 
to find $ (B) . 



* 2 <B> 



[l 



AS 

4r (B)M 



+ 



iKM 

2/ l M 2 -l 




= * 2 (C) 




AS iKM 

4r (C) M 2/M"-l 




(B) 

(C) 



‘ AS 
4r (B) M 






AS 

4r(C)M + 



iKM 
2/M 2 -1 



iKM 

2/h 2 -l 



AS 



AS 




AS l 
4 




( 3 . 20 ) 



*<c> [ k2 - Fn§T3i 2 ]^ + *< c > [ Ki - ] 4 

Then there are two equations [3.19 and 3.20] in two 

unknowns, and ^(B), which can be solved by using 

2 

Cramer's rule. The following abbreviations are defined as 



Crandall, S. H., Engineering Analysis, New York, 1956, 
p. 36. 
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A = 



A i “ 



A 2 = 



B = 



1 


AS 


iKM 


± 


4 r ( B ) M 


2/M2-1 AS 


1 + 


• AS 


iKM AC 

w=i A 


4r (A) M 


1 + 


AS 


- iKM AS 


4r (C)M 


2V^ 2 -1 



AS 



a=. + iKM ac - ' AS 2 , n 2 

4r(B)M 2/Mz-l 4 K ~ r 2 (B) M 2 

[- 



(3.21) 



C - * l( A, Al+ * 2 <A, L - + ^ 

AS + | *<A, [ K » - 



+ \ * (A) 



K 2 - 



r" (A)M' 



D = $2 (C)A 2+ * 1 (C) [- 4^ - + A|i (k* - F r T g T p-)] 



J t (C) [K- r ^]AS + I »(A) [k*. ^ 



2 “1 
(bTm : 



AS 



In terms of these abbreviations the two equations are 
A4> 1 (B) + B $ 2 (B) = C 

B$ 1 (B) + A $2 ( B ) = D (3.22) 

Then using Cramer 1 s rule 

2. ( B ) = CA - BD/A 2 - B 2 

$ 2 (B) = DA - CD/A 2 - B 2 - (3.23) 

was integrated to find 4>(B)along the right running char- 
acteristic and as a check $2 was also integrated along the 
left running characteristic using the trapezoidal rule 
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$ 1 (B) = $(A) + j j^>2 (A) + $ 2 (B) 1 AS 

$ 2 (B) = $(C) + | ^(C) + $ X (B)J AS 



(3.24) 



The difference 



was a measure of the 



$ 1 (B) - $ 2 (B) 

accuracy. If it becomes too large, the grid size must be 
made smaller. (B) and $ 2 (B) were averaged together for 
the best results. 



-MB) 




$ 1 (B) + $ 2 (B) 



(3.25) 



The above procedure will suffice for a general field 
point which is not subject to a boundary condition once the 
procedure has been initiated. 

On the initial Mach line it is known that $ = 0 , = 0, 

A$i = Of therefore from Eqs . 3.20 and 3.21 (See Fig. 5). 



a$ 2 (e) = a 2 $ 2 (d) 



2 

$2 (E) = $ 2 (D) A" 




(3.26) 
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At the initial point on the initial Mach line the panel 
boundary condition also applied, i.e., 



also from Eq. 3.14 

= 2 ( $ 1 “ *2) 

remembering that = 0 and combining the two equations 

# 2 = - — - iKZ 

2 M [ax J 

On the panel from the above boundary condition (See 
Fig. 6) 



$ 1 (L) - $ 2 (L) = - 1 |j| + iKZ J ( 3.2 

and from Eqs. 3.19 and 3.21 $^(L) + ^(L) was found on the 
right running characteristic and once again there exists a 
system of two equations in two unknowns. 





(3.29) 



A$ ^ (L) + B$ 2 (h) — C 



Using Cramer's rule 





(3.30) 



3 



Ibid. 
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FIGURE 6 

To write the program for a vibrating panel or the two- 
dimensional case, the radius of the cylinder r was set to 
infinity. Inspection of the basic equations programmed 
(Eq. 3.21) showed that r was always in the denominator, 
which causes these terms to equal zero. Therefore, terms 
containing r were neglected in the Fortran IV encoding of 
these equations. For the vibrating shell or three dimensional 
case the full equations as developed were Fortran IV encoded. 

Also, it should be noted that the problem was formulated 
for an arbitrary wall deflection Z (x) . However, in the 
numerical computations only the case of a sinsuoidal standing 
wave given by 



Z (x) = Z Q Sin (rruTx) (3.31) 

was considered. The axial wave-length 1 is hence equal to 
L/2m where m is the number of axial half-waves in a given 
cylinder length L. 
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C. TEST CASES 



The first test case considered a two-dimensional flat 
panel. If the cylinder radius is set at a very large value 
in relation to the length, then one has essentially a flat 
panel in supersonic flow. 




When MFREQ is set at various values, NFREQ and K set at zero, 
one then has supersonic flow over a wavy wall as developed 
in Johns Ref. 1 . 




Taking the wavy wall equations for $, <J> r and $ x from Johns 
Ref. 1 and normalizing results in 



$ = - TmTTY sin < x “ Jm 2 -1 rj 

$ r = 7T COS £ TT (X - j/S^-l rj^ 

= - /m 2 -” COS [* ( X " r ). 



( 3 . 31 ) 
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By using a large value of R in the three-dimensional case, 
it and the two-dimensional case results were then compared 
to the wavy wall results. These results were found to agree 
to the fifth significant figure. This case was also used ' 
to remove the majority of minor programming errors. 

To remove any other hidden programming errors it was 
decided to set each variable equal to a finite value and to 
step through the first three characteristic lines by hand. 
The values were then compared to those computed by the 
computer. In this manner the rest of the programming bugs 
were found and removed. 
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IV. COMPUTER PROGRAM ORGANIZATION AND OPERATION 



A. GENERAL ORGANIZATION 

The numerical methods discussed in the foregoing 
sections have been coded as a FORTRAN IV computer program 
for the IBM 360 system. 

The program consists of a main control section and four 
subsections or subroutines. An algorithm showing the 
program information flow and the interaction between the 
control section and the subroutines was listed in Figure 7. 
The control section reads in the input data, computes the 
characteristic grid, computes the initial condition, and 
computes several constants common to the subroutines. Then 
subroutine COMPRX is called to compute the coordinates 
of the next grid point to be calculated, see Figure 8 for 
a diagram of the characteristic grid used. Two "logical if 
statements then determine if the point in question is 
located on either the initial Mach line, the panel, or at 
a general field point. Control is then passed to that 
subroutine to compute the velocity potential and its 
derivatives along the characteristics at that point. The 
coefficient of pressure is computed along the panel in 
subroutine PANPT and after the program has stepped down to 
the last point on the last characteristic, a listing of the 
coefficient of pressure at each panel point is printed out. 
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COMPUTE 

CONSTANTS, 

SET 

COUNTERS ' 





CALL 

MACHLN 





■> 



> 




Yes 



RESET 
FLAGS & 
COUNTERS 



X 
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B. MAIN CONTROL PROGRAM OPERATION 



The main program dimensions the storage requirements 
in common block form; then it reads three input data cards. 

The first card has the date of the run in 3A4 format. The 
second card reads in the fineness of the grid, number of 
this run, circumferential mode number, and the axial mode 
number in 415 format. The third card reads Mach number, 
the reduced frequency number, and the cylinder radius in 
3F10.4 format. Input is then written to verify input data. 

The main program then computes the Mach angle, step 
sizes, and zeros the counters. The values for the velocity 
potential and its derivatives and the coordinates of each 
point are stored as a n x 2 matrix. As the program steps 
from characteristic- line to characteristic line, only those 
values for the line being computed and the line just 
computed are in storage. The initial values at the first 
point on the initial characteristic is computed using Eq. 

3.27 and written out to aid in the program check-out. The 
line counter, "Kount" , is then set at two and subroutine 
COMPXR is called to compute the values of x and r at the 
next point. The main program next computes various quantities 
dependent on x and r which are needed in the other subroutines 
The values for x, r, and dependent quantities are printed out 
as a program debug aid. Two "logical if" statements next 
determine if the point under consideration is on the initial 
Mach line, at a general field point, or on the panel. The 
appropriate subroutine is then called to solve for the . ' 
velocity potential and its derivatives along the 
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characteristics at that point. The next "logical if" 
statement determines if the program must switch lines 
or step along the present line. Two more "logical if" 
statements are used to determine if this solution is 
completed and if there is more data to be read and more 
solutions to be computed. If "numrun" is set at two, 
the program will look for more data; if "numrun" equals 
one, the program will terminate. Before it stops or 
reads more data the coefficient of pressure along the 
panel will be printed out for each panel point and an in- 
line graph, using the IBM scientific library program POLTP , 
will plot the real and imaginary parts of the coefficient 
of pressure. 

C. SUBROUTINE COMPXR OPERATION 

Information is passed between the calling program and 
the subroutines by means of the common storage. The first 
step in COMPXR is to determine if the point to be computed 
is on the initial Mach line or not. The coordinates of the 
last point computed are used and the step distance in the 
x and r directions either added to or subtracted from them 
to find the new x and r coordinates. The addition or sub- 
traction depends upon whether we step up to a new character- 
istic line or down the present characteristic line. 
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D. SUBROUTINE MACHLN OPERATION 

In subroutine MACHLN the boundary conditions of the 
equations developed in the foregoing section are enforced 
such that the velocity potential and its derivatives along 
the right running characteristic are set to zero. Equations 
3.26 were coded in FORTRAN to solve for a t the point in 

question. The value of $ 2 and the variables used to solve 
for $2 are printed out for the first 12 characteristic lines 
to aid in the program checkout and then control is returned 
to the main program. 

E. OPERATION OF SUBROUTINE PANPT 

In subroutine PANPT the boundary condition along the 
panel is enforced. The velocity potential and its derivatives 
along the characteristics are computed at the point in 
question, using Eq. 3.29. The real and imaginary parts of 
the coefficient of pressure are stored for future read out. 
Control is then returned to the main program. 

F. OPERATION OF SUBROUTINE GENFPT 

GENFPT computes the velocity potential and its first 
derivatives along the right and left running characteristics. 
It draws the necessary information from common storage and 
uses the Eqs. 3.21, 3.22, and 3.23, which were developed in 
the foregoing chapter. To aid in program check out the 
variables A, A 1 , A 2 , B, C, and D are printed out for the 
first 12 characteristic lines. Control is then returned to 
the main program. 
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V. RESULTS AND DISCUSSION 



The case of a vibrating two-dimensional panel was 
chosen as the first test of the accuracy of the character- 
istics method. These results were compared to the results 
obtained by evaluating the method-of-singularities integrals 
developed in Ref. 4 for flutter of a two-dimensional panel 
with one surface exposed to supersonic flow. These integrals 
were evaluated by applying a 12-point Gaussian Quadrature. 

The results of these evaluations were available in unpublished 
notes by M. F. Platzer. 

The comparisons for M = v'§~ at a reduced frequency of 0.5 

and an axial mode number of one are shown in Figures 8-11. 
Figures 8 and 9 are plots of the real and imaginary parts of 
the coefficient of pressure obtained by the method of singu- 
larities. Figures 10 and 11 are plots of the real and 
imaginary parts of the Cp for the characteristics method. 

When these tests were first run it was found that the real 
parts of the Cp plots were identical, as shown in Figures 8 
and 10. The imaginary parts of Cp, shown in Figures 9 and 11, 
did not agree. In rechecking the formulation of the equations 
it was found that the values of $ , and <J> 2 had not been 
averaged as explained in Chapter II, Section B, where Eqs. 

3.19 and 3.20 were formulated. When these average values 
were incorporated into the equations, the computer program 
rewritten and the program rerun, then the two plots were 

► " 

found to be identical for both the real and imaginary parts 



of the Cp. 
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Figures 12-14 are plots for a reduced frequency of two, 



an axial mode of four, and the same Mach number. These 
comparisons of the real and imaginary parts of the Cp were 
also found to be identical. 

In attempting to verify the results of the vibrating 
shell program, Ref. 5 was used. Reference 5 includes a 
study of supersonic potential flow past an infinitely long 
vibrating cylinder. It was not known what effect the dif- 
ference in an infinitely long cylinder vice a finite length 
cylinder would have on the outcome, but these were the only 
published results found that could be used for a comparison 
of the vibrating shell program. 

Two parameters were used to compare the results of the 
two studies: the maximum amplitude, A^, of the sinsuoidal 

Cp wave, and the phase angle, of this wave. Equation 14 

in Ref. 5 is used to equate Cp to the value of A-^ listed in 
Figure 2a of Ref. 5. 



— Z o pU 2 mirR 

p(x,R,e) = R~ -jf- “l” Cos n6 A x Cos 




p is defined as p = p - p ro where 



P " 



P CX) 



(5.1) 



C P - % pU 



43 



and the terms, Z L, and Cos(x) were all taken as equal 
to one then 

C n = -r — 2^ = 2. mirA 

P h pU" M 1 

A 1 = C p (5.2) 

Three different values of n were computed and the coefficients 
of pressure were listed in Figure 17 as a test case for this 
comparison. For the case n = 0, (See Figure 2a of Ref. 5) 
we find that for the ratio of axial wave length/cylinder 
radius £/R = 3.33, A^ equals 1.05. Using Eq. 5.2 to determine 
A^ for a maximum Cp of 5.63 in Figure 17, A^ for this study 
would equal 1.044. Using this same procedure for n = 4 and 8 

n = 4 A 1 (ref 5) = 1.30 A = 1.445 

n = 8 A^(ref 5) = 1.80 A^ = 1.85. 

To compare the phase angle ^ , Figure 2b in Ref. 5 was 
entered for %/ R = 3.33 and the three values of n taken 
previously. This shift, was calculated from results 

(Fig. 17) by determining the zero crossing point. 

n = 0 59.8 

n = 4 * 57.7 

n = 8 48.9 

Converting these zero crossing points to angular shifts, 
it was found that 
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n = 


0 


=• 4.8 


*1 


(ref 


5) = 5 


n = 


4 


i> 1 = 13.5 


*1 


(ref 


5) = 10 


n = 


8 


= 55.4 


*1 


(ref 


5) = 60 



which compare favorably with Anderson's values listed in 
the third column. 



45 



m =J, K =0.5 
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METHOD OF SINGULARITIES 



so 
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METHOD OF SINGULARITIES 
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METHOD OF CHARACTERISTICS 
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METHOD OF CHARACTERISTICS 
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METHOD OF SINGULARITIES 
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METHOD OF SINGULARITIES 
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METHOD OF CHARACTERISTICS 
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METHOD OF CHARACTERISTICS 



COEFFICIENT OF PRESSURE COMPUTER OUTPUT LISTING 



M 

where 



FOR CASE OF VIBRATING SHELL 



= 3.5 
L is 


, m = 3 , K = 0 , r 
divided into 120 


= 0.20 and n = 
grid points. 


0, 1, 3, 


Grid 


Point 


n = 0 


n = 4 


n = 8 


1 . 


5.6198 


5.6198 


5.6198 


2. 


5.5677 


5.5538 


5.5642 


3. 


5.4813 


5.4266 


5.4676 


4. 


5.3613 


5.2396 


5.3307 


5. 


5.2083 


4.9950 


5.1546 


6 • 


5.0233 


4.6954 


4.9404 


7. 


4.8075 


4.3440 


4.6897 


8. 


4.5621 


3.9444 


4.4041 


9. 


4.2887 


3.5007 


4.0857 


10. 


3.9890 


3.0172 


3.7366 


11. 


3.6647 


2.4988 


3.3591 


12. 


3.3180 


1.9504 


2.9559 


13. 


2.9509 


1.3775 


2.5295 


14. 


2.5657 


0.7855 


2.0828 


15. 


2.1648 


0.1801 


1.6188 


16. 


1.7507 


-0.4329 


1.1407 


17. 


1.3258 


-1.0477 


0.6515 


18. 


0.8929 


-1.6586 


0.1546 


19. 


0.4546 


-2.2597 


-0.3467 


20. 


0.0136 


-2.8453 


-0.8492 


21. 


-0.4273 


-3.4100 


-1.3495 


22. 


— 0 . 8656 


-3.9435 


-1.8442 


23. 


-1.2985 


-4.4557 


-2.3301 


24. 


-1.7232 


-4.9267 


-2.8040 


25. 


-2.1373 


-5.3572 


-3.2626 


26. 


-2.5380 


-5.7431 


-3.7028 


27. 


-2.9230 


-6.0806 


-4.1218 


28. 


-3.2899 


-6.3665 


-4.5167 


29. 


-3.6365 


‘ -6.5980 


-4.8848 


30. 


-3.9605 


-6.7726 


-5.2237 


31. 


-4.2600 


-6.8886 


-5.5309 


32. 


-4.5331 


-6.9446 


-5.8044 


33. 


-4.7782 


-6.9398 


-6.0423 


34 . 


-4.9938 


-6.8738 


-6.2428 


35. 


-5.1785 


-6.7468 


-6.4046 


36. 


-5.3312 


-6.5595 


-6.5263 


37. 


-5.4509 


-6.3132 


-6.6071 


38. 


-5.5369 


-6.0095 


-6.6462 


39. 


-5.5887 


-5.6507 


-6.6432 
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Grid 



Point 


n = 0 


n = 4 


3 

t! 

00 


40. 


-5.6059 


-6.5979 


-5.2394 


41. 


-5.5885 


-6.5104 


-4.7787 


42. 


-5.5366 


-6.3810 


-4.2721 


43. 


-5.4504 


-6.2105 


- 3.7235 


44. 


-5.3305 


-5.9996 


-3.1371 


45. 


-5.1777 


-5.7495 


-2.5175 


46. 


-4.9929 


-5.4616 


-1.8695 


47. 


-4.7773 


-5.1375 


-1.1981 


48. 


-4.5320 


-4.7791 


-0.5086 


49. 


-4.2588 


-4.3884 


-0.1934 


50. 


-3.9592 


-3.9678 


-0.9027 


51. 


-3.6351 


-3.5196 


-1.6137 


52. 


-3.2886 


-3.0466 


-2.3205 


53. 


-2.9216 


-2.5516 


-3.0178 


54. 


-2.5366 


-2.0375 


-3.7000 


55. 


-2.1359 


-1.5074 


-4.3616 


56. 


-1.7219 


-0.9644 


-4.9975 


57. 


-1.2972 


-0.4119 


-5.6025 


58. 


-0.8645 


0.1467 


-6.1719 


59. 


-0.4263 


0.7082 


-6.7010 


60. 


0.0145 


1.2692 


-7.1856 


61. 


0.4554 


1.8261 


7.6218 


62. 


0.8935 


2.3757 


8.0060 


63. 


1.3262 


2.9145 


8.3350 


64. 


1.7508 


3.4393 


8.6062 


65. 


2.1647 


3.9467 


8.8171 


6 6 . 


2.5653 


4.4338 


8.9659 


67. 


2.9502 


4.8975 


9.0512 


68. 


3.3170 


5.3350 


9.0720 


69. 


3.6633 


5.7435 


9.0280 


70. 


3.9872 


6.1205 


8.9191 


71. 


4.2866 


6.4636 


8.7458 


72. 


4.5596 


6.7708 


8.5092 


73. 


4.8046 


7.0400 


8.2107 


74. 


5.0200 


7.2697 


7.8522 


75. 


5.2046 


7.4582 


7.4361 


76. 


5.3571 


7.6045 


6.9652 


77. 


5.4767 


7.7075 


6.4427 


78. 


5.5626 


7.7666 


5.8721 


79. 


5.6142 


7.7813 


5.2573 


80. 


5.6313 


7.7514 


4.6027 


81. 


5.6138 


' 7.6771 


3.9126 


82. 


5.5617 


7.5588 


3.1920 


83. 


5.4755 


7.3971 


2.4459 


84. 


5.3555 


7.1929 


1.6793 


85. 


5.2025 


6.9474 


0.8978 
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Grid 



Point 


3 

II 

o 


n = 4 


3 

II 

00 


86. 


5.0176 


6.6619 


0.1068 


87. 


4.8018 


6.3383 


-0.6882 


88. 


4.5565 


5.9783 


-1.4817 


89. 


4.2831 


5.5841 


-2.2681 


90. 


3.9834 


5.1580 


-3.0419 


91. 


3.6592 


4.7026 


-3.7977 


92. 


3.3125 


4.2205 


-4.5302 


93. 


2.9455 


3.7146 


-5.2342 


94. 


2.5604 


3.1879 


-5.9048 


95. 


2.1595 


2.6436 


-6.5374 


96. 


1.7454 


2.0848 


-7.1275 


97. 


1.3206 


1.5150 


-7.6710 


98. 


0.8878 


0.9376 


-8.1640 


99. 


0.4495 


0.3559 


-8.6031 


100. 


0.0085 


-0.2264 


-8.9853 


101. 


-0.4324 


-0.8061 


-9.3078 


102. 


-0.8706 


-1.3795 


-9.5684 


103. 


-1.3034 


-1.9434 


- 9.7652 


104. 


- 1.7282 


-2.4942 


-9.8968 


105. 


-2.1421 


-3.0288 


-9.9623 


106. 


-2.5429 


-3.5439 


-9.9611 


107. 


-2.9278 


-4.0364 


-9.8932 


108. 


-3.2947 


-4.5035 


-9.7591 


109. 


-3.6412 


-4.9423 


-9.5595 


110. 


-3.9652 


-5.3503 


-9.2957 


111. 


-4.2646 


- 5.7250 


-8.9696 


112. 


-4.5378 


-6.0642 


-8.5833 


113. 


-4.7828 


-6.3660 


-8.1393 


114 . 


-4.9984 


-6.6285 


-7.6406 


115. 


-5.1830 


-6.8502 


-7.0906 


116. 


-5.3356 


- 7.0299 


-6.4928 


117. 


-5.4553 


-7.1665 


-5.8513 


118. 


-5.5413 


- 7.2593 


-5 ,,17 03 


119. 


-5.5931 


-7.3078 


-4.4545 


120. 


-5.6103 


-7.3117 


-3.7084 
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VI. SUMMARY AND CONCLUSIONS 



A linearized method of characteristics procedure was 
developed to compute supersonic flow past vibrating panels 
and shells. Good agreement was obtained when comparing 
with previous results by Nelson and Cunningham {Ref. 4] 
and Anderson [Ref. |] , who used quite different approaches 
to this problem. 

However, more work needs to be done to further verify 
the cylindrical shell computations. In particular, the 
computations of the generalized aerodynamic forces presented 
by Dowell and Widnall [kef. (T) should be verified. For this 
purpose the limiting case of vanishing shell radius (slender 
body theory) must be studied. This requires a reformulation 
of the flow boundary condition. 

Also, it should be noted that the characteristics method 
allows the incorporation of quite arbitrary vibration modes 
(not restricted to sinsuoidal modes) and hence should prove 
to be a quite versatile tool for aeroelastic analyses in the 
supersonic flow regime. 
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APPENDIX A: VIBRATING PANEL COMPUTER P ROGRAM 
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PRINT NAME OF PROGRAM AND DATE OF RUN 
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LINEARIZED METHOD OF CHARACTERISTICS 
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COEFFICENT 


OF PRESSURE 




1.0 


0.0 


0.0 


2.0 


24.55586270 


0.00408154 


3.0 


22. 85066452 


0.01351091 


4.0 


20.09254127 


0.03334298 


5.0 


16. 40344123 


0. 06806011 


6.0 


11.94654108 


0. 12131901 


7.0 


6. 9190 92 87 


0.19573567 


8.0 


1. 54377602 


0.29271826 


0.0 


-3.94106736 


0.41235615 


10.0 


-9. 29188676 


0.55337020 


11.0 


-14. 27060977 


0. 71 312712 


12.0 


-18.65509318 


0.88771781 


13.0 


-22. 24888160 


1. 07209644 


14.0 


-24.88984639 


1.26027485 


15.0 


-26.45732668 


1.44556368 


16.0 


-26. 87745992 


1.62085036 


17.0 


-26.12646822 


1.77890147 


18.0 


-24. 23175624 


1.91267680 


19.0 


-21.27077195 


2. 01564102 


20.0 


-17.36767944 


2.08205951 


21.0 


- 12. 68798847 


2. 10726521 


22.0 


-7.43137495 


2.08788457 


23.0 


-1. 823C0542 


2.02201259 


24.0 


3. 89625584 


1. 90932880 


25.0 


9.48032952 


1.75114890 


26.0 


14.68863637 


1.55040919 


27.0 


19.29657652 


1. 31158411 


28.0 


23.10531564 


1.04053992 


29.0 


25. 95045C47 


0. 74433023 


30.0 


27.7091 7520 


0. 43094193 


31.0 


28.30563578 


0.10900179 


32.0 


27. 71423823 


-0. 21254394 


33.0 


25.96076634 


-0.52476371 


34.0 


23. 12126002 


-0.81907450 


35.0 


19.31870316 


-1.08756155 


36.0 


14.71766553 


-1.32327050 


37.0 


9. 51713284 


-1. 52045979 


38.0 


3.94183787 


-1. 67480302 


39.0 


-1. 76752874 


-1.78353270 
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40.0 


- 7. 36480060 


- 1 . 84551983 


41.0 


- 12.60 905291 


- 1. 86128623 


42.0 


- 17.27508720 


- 1.83294939 


43.0 


- 21 . 16322567 


- 1 . 76410272 


44.0 


- 24.10798772 


- 1.65963661 


45.0 


- 25.98527008 


- 1 . 52550839 


46.0 


- 26.71771716 


- 1 . 36847141 


47.0 


- 26.27804790 


- 1.19577514 


48.0 


- 24.69019390 


- 1. 01484935 


49.0 


- 22.02820014 


- 0. 83298627 


50.0 


- 18.41293660 


- 0.65703434 


51.0 


- 14.00676536 


- 0.49311690 


52.0 


- 9.00639664 


- 0.34638770 


53.0 


- 3. 63424668 


- 0.22083373 


54.0 


1.87132446 


- 0. 11913365 


55.0 


7.26646797 


- 0.04257745 


56.0 


12 . 31253691 


0.00894953 


57.0 


16.78654483 


0.03691793 


58.0 


20.49093494 


0.04404391 


59.0 


23. 26223271 


0.03412212 


60.0 


24.97820347 


0.01181008 



83 




84 




85 



LINEARIZED METHOD OF CHARACTERISTICS 



APPENDIX D: VIBRATING SHELL SAMPLE COMPUTER OUTPUT 
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IHAVEP FOR RUN =1 KQUNT = 4 

PX = -2.83598506 0.0 PR = 9. 51218306 
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